# input: fun,gfun,Hfun,x0
# output x, val , k
#mark = 100
#k = 0
#epsilon = 1e-5
#while k < mark :
#    gk= feval(gfun,x0)
#    Hk = feval(Hess,x0)
#    dk = -inv(Hk)\gk
#    if ( norm(gk)< epsilon):
#        break
#    x0 = x0 + dk
#    k=k+1
#
#x = x0
#val = feval(func,x)


def newton(fun,gfun,Hfun,x0,k,eps,maxk=100):
    gk= gfun(x0)
    Hk = Hfun(x0)
    if np.lingla.norm(gk) < eps or k >= maxk:
        return x0, fun(x0), k
    dk = -np.lingla.inv(Hk)\gk 
    return(fun,gfun,Hfun,x0+dk,k+1,eps,maxk)
